function f = fitness(Phen,expPara,modPara)

NIND = size(Phen,1);
f = zeros(NIND,2);

for i = 1:NIND;
    clf;

    for j=1:2
    covfunc = expPara{j,1};
    loghyper1 = expPara{j,2};
    loghyper2 = expPara{j,3};
    tmp = expPara{j,4};

    x  = tmp(:,1);
    y1 = tmp(:,2);
    y2 = tmp(:,3);

    sample = modPara{j,1};
    column = modPara{j,2};
    elution = modPara{j,3};
    ocfe = modPara{j,4};
        
    isotherm = Phen(i,1:7);
    if (j==2);
        isotherm(7) = Phen(i,8);
    end
    load VS_DL.mat
    [DL1] = gpr(loghyper_BSA, covfunc, VS_BSA*v, DL_BSA*v, isotherm(7)*v);
    [DL2] = gpr(loghyper_MYO, covfunc, VS_MYO*v, DL_MYO*v, isotherm(7)*v);
    
    isotherm(8) = DL1/v;
    isotherm(9) = DL2/v;
    
    [ts, C, C1, C2] = ED_GRsolver(sample, column, elution, isotherm, ocfe);
    
    %fprintf('DL: %f, DL1: %f, DL2: %f\n', DL, DL1, DL2);
    
    tn = 1000;
    cv = 10^6;
    tstar = ts/tn;    
    cv1 = 7.5639e-6*cv;
    cv2 = 5.86407e-6*cv;
    C1 = C1/sample(2)*cv1;
    C2 = C2/sample(3)*cv2;

    [mu S2] = gpr(loghyper1, covfunc, x, y1, tstar);
    %mu_M = mu;
    %S2_M = S2;

    clf;
    hold on
    f1 = Objfun(C1, mu, S2);
    tstr = sprintf('C1 f:%f',f1);
    subplot(2,1,1);
    drawFig(x,y1,mu,C1,S2,tstar,tstr);

    [mu S2] = gpr(loghyper2, covfunc, x, y2, tstar);
    %mu_B = mu;
    %S2_B = S2;
         
    f2 = Objfun(C2, mu, S2);
    tstr = sprintf('C2 f:%f',f2);
    subplot(2,1,2);
    drawFig(x,y2, mu,C2,S2,tstar,tstr);
    hold off
    %     C = [C1; C2];
    %     mu = [mu_M; mu_B];
    %     S2 = [S2_M; S2_B];
    %     f(i) = Objfun(C, mu, S2);
    
    if f1==1000 && f2==1000;
        f(i,j) = 1000;
    elseif f1==1000 || f2==1000;
        f(i,j) = 999;
    else
        f(i,j) = max(f1,f2);
    end
    fprintf('Num: %d, Daten: %d, f1: %f, f2: %f, f: %f\n',i,j,f1,f2,f(i,j));
    
   % pause;
    
    %f(i) = f1;
    %fprintf('Num: %d, f: %f\n',i,f(i));
    %fprintf('Lambda: %f, K1: %f, K2: %f, sigma1: %f, sigma2: %f, nu1: %f, nu2: %f', isotherm);
    
    end
end

f = f(:,1)+f(:,2);